 ! #############################################################################
 ! Form the antisymmetrised spinorbital basis ERIs: <pq||rs> = <pq|rs> - <pq|sr>
 ! where <pq|rs> = (pr|qs) * delta(sigma_p,sigma_r) * delta(sigma_q,sigma_s)
 ! if sigma_p and sigma_r are the same ,delta value is 1, different value is 0
 !该子程序的双电子积分用的化学符号表示方法
 ! #############################################################################
    subroutine MO2SOERI(NORB,ERIMO,ERISPIN)
        implicit none
        integer :: NORB,nspin
        integer :: ip,iq,ir,is
        integer :: ipspa,iqspa,irspa,isspa
        integer :: modip2,modiq2,modir2,modis2
        real*8  :: V1,V2
        real*8  :: krodelta 
        real*8  :: kdeltpr,kdeltqs,kdeltps,kdeltrq
        real*8 :: ERIMO(NORB,NORB,NORB,NORB)
        real*8 :: ERISPIN(2*NORB,2*NORB,2*NORB,2*NORB)
!!!!!!!!     <pq||rs>
!!!!!!!!    i,j,k,l is MO ,  p,q,r,s is  spin-orbital.
        write(21,*)'==========================================='
        write(21,*)'   *** Start MO to spin-orbital ERI ***' 
        write(21,*)'   ***  MO to spin-orbital  *** '       
        !write(*,*)ERIMO(1,1,1,1)
        ERISPIN=0d0
        nspin=2*NORB
        do ip=1,nspin
            ipspa=ceiling(ip/2d0)
            modip2=mod(ip,2)
            do iq=1,nspin  
            iqspa=ceiling(iq/2d0)
            modiq2=mod(iq,2) 
            !kdeltpq= krodelta(modip2,modiq2)
                do ir=1,nspin    
                    irspa=ceiling(ir/2d0)
                    modir2=mod(ir,2)
                    do is=1,nspin
                        isspa=ceiling(is/2d0)
                        modis2=mod(is,2)
                        kdeltpr= krodelta(modip2,modir2)
                        kdeltqs= krodelta(modiq2,modis2)
                        kdeltps= krodelta(modis2,modip2)
                        kdeltrq= krodelta(modir2,modiq2)
!!!!                   
                        v1=ERIMO(ipspa,irspa,iqspa,isspa)*kdeltpr*kdeltqs
                        v2=ERIMO(ipspa,isspa,irspa,iqspa)*kdeltps*kdeltrq
                        ERISPIN(ip,iq,ir,is)=v1-v2    !<>
 !                       write(25,*)ip,iq,ir,is,ERISPIN(ip,iq,ir,is)
                    end do
                end do 
            end do 
        end do 
        write(21,*)'   ***  MO to spin-orbital Done ***' 
        write(21,*)'==========================================='   
      
               
    end subroutine    
!!!!===============================================================!!!!!!!!!!
     real*8 function krodelta(i,j)
!!!! if i=j  ,delta =1
!!!! if i/=j ,delta =0 
            implicit none
            integer :: i,j
            if(i==j)then
               krodelta=1
            else
               krodelta=0
            end if
            return
        end function    krodelta

  !该子程序的双电子积分用的化学符号表示方法
    subroutine MO2SOERI1(NORB,ERIMO,ERISPIN)
        implicit none
        integer :: NORB,nspin
        integer :: ip,iq,ir,is
        integer :: ipspa,iqspa,irspa,isspa
        integer :: modip2,modiq2,modir2,modis2
        real*8  :: V1,V2
        real*8  :: krodelta
        real*8  :: kdeltpq,kdeltrs,kdeltps,kdeltrq
        real*8 :: ERIMO(NORB,NORB,NORB,NORB)
        real*8 :: ERISPIN(2*NORB,2*NORB,2*NORB,2*NORB)
!!!!!!!!   (ij||kl)  -> (pq||rs)
!!!!!!!!    i,j,k,l is MO ,  p,q,r,s is  spin-orbital.
        ERISPIN=0d0
        nspin=2*NORB
        do ip=1,nspin
            ipspa=ceiling(ip/2d0)
            modip2=mod(ip,2)
            do iq=1,nspin
            iqspa=ceiling(iq/2d0)
            modiq2=mod(iq,2)
            kdeltpq= krodelta(modip2,modiq2)
                do ir=1,nspin
                    irspa=ceiling(ir/2d0)
                    modir2=mod(ir,2)
                    do is=1,nspin
                        isspa=ceiling(is/2d0)
                        modis2=mod(is,2)
                        kdeltrs= krodelta(modis2,modir2)
                        kdeltps= krodelta(modis2,modip2)
                        kdeltrq= krodelta(modir2,modiq2)
!!!!                    
                        v1=ERIMO(ipspa,iqspa,irspa,isspa)*kdeltpq*kdeltrs
                        v2=ERIMO(ipspa,isspa,irspa,iqspa)*kdeltps*kdeltrq
                        ERISPIN(ip,iq,ir,is)=v1-v2
!                        write(26,*)ip,iq,ir,is,ERISPIN(ip,ir,iq,is)
                    end do
                end do
            end do
        end do
    end subroutine

